home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Cream of the Crop 26
/
Cream of the Crop 26.iso
/
os2
/
octa209s.zip
/
octave-2.09
/
libcruft
/
minpack
/
hybrj.f
< prev
next >
Wrap
Text File
|
1996-07-19
|
14KB
|
442 lines
SUBROUTINE HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,DIAG,MODE,
* FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF,WA1,WA2,
* WA3,WA4)
INTEGER N,LDFJAC,MAXFEV,MODE,NPRINT,INFO,NFEV,NJEV,LR
DOUBLE PRECISION XTOL,FACTOR
DOUBLE PRECISION X(N),FVEC(N),FJAC(LDFJAC,N),DIAG(N),R(LR),
* QTF(N),WA1(N),WA2(N),WA3(N),WA4(N)
EXTERNAL FCN
C **********
C
C SUBROUTINE HYBRJ
C
C THE PURPOSE OF HYBRJ IS TO FIND A ZERO OF A SYSTEM OF
C N NONLINEAR FUNCTIONS IN N VARIABLES BY A MODIFICATION
C OF THE POWELL HYBRID METHOD. THE USER MUST PROVIDE A
C SUBROUTINE WHICH CALCULATES THE FUNCTIONS AND THE JACOBIAN.
C
C THE SUBROUTINE STATEMENT IS
C
C SUBROUTINE HYBRJ(FCN,N,X,FVEC,FJAC,LDFJAC,XTOL,MAXFEV,DIAG,
C MODE,FACTOR,NPRINT,INFO,NFEV,NJEV,R,LR,QTF,
C WA1,WA2,WA3,WA4)
C
C WHERE
C
C FCN IS THE NAME OF THE USER-SUPPLIED SUBROUTINE WHICH
C CALCULATES THE FUNCTIONS AND THE JACOBIAN. FCN MUST
C BE DECLARED IN AN EXTERNAL STATEMENT IN THE USER
C CALLING PROGRAM, AND SHOULD BE WRITTEN AS FOLLOWS.
C
C SUBROUTINE FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
C INTEGER N,LDFJAC,IFLAG
C DOUBLE PRECISION X(N),FVEC(N),FJAC(LDFJAC,N)
C ----------
C IF IFLAG = 1 CALCULATE THE FUNCTIONS AT X AND
C RETURN THIS VECTOR IN FVEC. DO NOT ALTER FJAC.
C IF IFLAG = 2 CALCULATE THE JACOBIAN AT X AND
C RETURN THIS MATRIX IN FJAC. DO NOT ALTER FVEC.
C ---------
C RETURN
C END
C
C THE VALUE OF IFLAG SHOULD NOT BE CHANGED BY FCN UNLESS
C THE USER WANTS TO TERMINATE EXECUTION OF HYBRJ.
C IN THIS CASE SET IFLAG TO A NEGATIVE INTEGER.
C
C N IS A POSITIVE INTEGER INPUT VARIABLE SET TO THE NUMBER
C OF FUNCTIONS AND VARIABLES.
C
C X IS AN ARRAY OF LENGTH N. ON INPUT X MUST CONTAIN
C AN INITIAL ESTIMATE OF THE SOLUTION VECTOR. ON OUTPUT X
C CONTAINS THE FINAL ESTIMATE OF THE SOLUTION VECTOR.
C
C FVEC IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS
C THE FUNCTIONS EVALUATED AT THE OUTPUT X.
C
C FJAC IS AN OUTPUT N BY N ARRAY WHICH CONTAINS THE
C ORTHOGONAL MATRIX Q PRODUCED BY THE QR FACTORIZATION
C OF THE FINAL APPROXIMATE JACOBIAN.
C
C LDFJAC IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN N
C WHICH SPECIFIES THE LEADING DIMENSION OF THE ARRAY FJAC.
C
C XTOL IS A NONNEGATIVE INPUT VARIABLE. TERMINATION
C OCCURS WHEN THE RELATIVE ERROR BETWEEN TWO CONSECUTIVE
C ITERATES IS AT MOST XTOL.
C
C MAXFEV IS A POSITIVE INTEGER INPUT VARIABLE. TERMINATION
C OCCURS WHEN THE NUMBER OF CALLS TO FCN WITH IFLAG = 1
C HAS REACHED MAXFEV.
C
C DIAG IS AN ARRAY OF LENGTH N. IF MODE = 1 (SEE
C BELOW), DIAG IS INTERNALLY SET. IF MODE = 2, DIAG
C MUST CONTAIN POSITIVE ENTRIES THAT SERVE AS IMPLICIT
C (MULTIPLICATIVE) SCALE FACTORS FOR THE VARIABLES.
C
C MODE IS AN INTEGER INPUT VARIABLE. IF MODE = 1, THE
C VARIABLES WILL BE SCALED INTERNALLY. IF MODE = 2,
C THE SCALING IS SPECIFIED BY THE INPUT DIAG. OTHER
C VALUES OF MODE ARE EQUIVALENT TO MODE = 1.
C
C FACTOR IS A POSITIVE INPUT VARIABLE USED IN DETERMINING THE
C INITIAL STEP BOUND. THIS BOUND IS SET TO THE PRODUCT OF
C FACTOR AND THE EUCLIDEAN NORM OF DIAG*X IF NONZERO, OR ELSE
C TO FACTOR ITSELF. IN MOST CASES FACTOR SHOULD LIE IN THE
C INTERVAL (.1,100.). 100. IS A GENERALLY RECOMMENDED VALUE.
C
C NPRINT IS AN INTEGER INPUT VARIABLE THAT ENABLES CONTROLLED
C PRINTING OF ITERATES IF IT IS POSITIVE. IN THIS CASE,
C FCN IS CALLED WITH IFLAG = 0 AT THE BEGINNING OF THE FIRST
C ITERATION AND EVERY NPRINT ITERATIONS THEREAFTER AND
C IMMEDIATELY PRIOR TO RETURN, WITH X AND FVEC AVAILABLE
C FOR PRINTING. FVEC AND FJAC SHOULD NOT BE ALTERED.
C IF NPRINT IS NOT POSITIVE, NO SPECIAL CALLS OF FCN
C WITH IFLAG = 0 ARE MADE.
C
C INFO IS AN INTEGER OUTPUT VARIABLE. IF THE USER HAS
C TERMINATED EXECUTION, INFO IS SET TO THE (NEGATIVE)
C VALUE OF IFLAG. SEE DESCRIPTION OF FCN. OTHERWISE,
C INFO IS SET AS FOLLOWS.
C
C INFO = 0 IMPROPER INPUT PARAMETERS.
C
C INFO = 1 RELATIVE ERROR BETWEEN TWO CONSECUTIVE ITERATES
C IS AT MOST TOL.
C
C INFO = 2 NUMBER OF CALLS TO FCN WITH IFLAG = 1 HAS
C REACHED MAXFEV.
C
C INFO = 3 XTOL IS TOO SMALL. NO FURTHER IMPROVEMENT IN
C THE APPROXIMATE SOLUTION X IS POSSIBLE.
C
C INFO = 4 ITERATION IS NOT MAKING GOOD PROGRESS, AS
C MEASURED BY THE IMPROVEMENT FROM THE LAST
C FIVE JACOBIAN EVALUATIONS.
C
C INFO = 5 ITERATION IS NOT MAKING GOOD PROGRESS, AS
C MEASURED BY THE IMPROVEMENT FROM THE LAST
C TEN ITERATIONS.
C
C NFEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
C CALLS TO FCN WITH IFLAG = 1.
C
C NJEV IS AN INTEGER OUTPUT VARIABLE SET TO THE NUMBER OF
C CALLS TO FCN WITH IFLAG = 2.
C
C R IS AN OUTPUT ARRAY OF LENGTH LR WHICH CONTAINS THE
C UPPER TRIANGULAR MATRIX PRODUCED BY THE QR FACTORIZATION
C OF THE FINAL APPROXIMATE JACOBIAN, STORED ROWWISE.
C
C LR IS A POSITIVE INTEGER INPUT VARIABLE NOT LESS THAN
C (N*(N+1))/2.
C
C QTF IS AN OUTPUT ARRAY OF LENGTH N WHICH CONTAINS
C THE VECTOR (Q TRANSPOSE)*FVEC.
C
C WA1, WA2, WA3, AND WA4 ARE WORK ARRAYS OF LENGTH N.
C
C SUBPROGRAMS CALLED
C
C USER-SUPPLIED ...... FCN
C
C MINPACK-SUPPLIED ... DOGLEG,DPMPAR,ENORM,
C QFORM,QRFAC,R1MPYQ,R1UPDT
C
C FORTRAN-SUPPLIED ... DABS,DMAX1,DMIN1,MOD
C
C MINPACK. VERSION OF SEPTEMBER 1979.
C BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
C
C **********
INTEGER I,IFLAG,ITER,J,JM1,L,NCFAIL,NCSUC,NSLOW1,NSLOW2
INTEGER IWA(1)
LOGICAL JEVAL,SING
DOUBLE PRECISION ACTRED,DELTA,EPSMCH,FNORM,FNORM1,ONE,PNORM,
* PRERED,P1,P5,P001,P0001,RATIO,SUM,TEMP,XNORM,
* ZERO
DOUBLE PRECISION DPMPAR,ENORM
DATA ONE,P1,P5,P001,P0001,ZERO
* /1.0D0,1.0D-1,5.0D-1,1.0D-3,1.0D-4,0.0D0/
C
C EPSMCH IS THE MACHINE PRECISION.
C
EPSMCH = DPMPAR(1)
C
INFO = 0
IFLAG = 0
NFEV = 0
NJEV = 0
C
C CHECK THE INPUT PARAMETERS FOR ERRORS.
C
IF (N .LE. 0 .OR. LDFJAC .LT. N .OR. XTOL .LT. ZERO
* .OR. MAXFEV .LE. 0 .OR. FACTOR .LE. ZERO
* .OR. LR .LT. (N*(N + 1))/2) GO TO 300
IF (MODE .NE. 2) GO TO 20
DO 10 J = 1, N
IF (DIAG(J) .LE. ZERO) GO TO 300
10 CONTINUE
20 CONTINUE
C
C EVALUATE THE FUNCTION AT THE STARTING POINT
C AND CALCULATE ITS NORM.
C
IFLAG = 1
CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
NFEV = 1
IF (IFLAG .LT. 0) GO TO 300
FNORM = ENORM(N,FVEC)
C
C INITIALIZE ITERATION COUNTER AND MONITORS.
C
ITER = 1
NCSUC = 0
NCFAIL = 0
NSLOW1 = 0
NSLOW2 = 0
C
C BEGINNING OF THE OUTER LOOP.
C
30 CONTINUE
JEVAL = .TRUE.
C
C CALCULATE THE JACOBIAN MATRIX.
C
IFLAG = 2
CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
NJEV = NJEV + 1
IF (IFLAG .LT. 0) GO TO 300
C
C COMPUTE THE QR FACTORIZATION OF THE JACOBIAN.
C
CALL QRFAC(N,N,FJAC,LDFJAC,.FALSE.,IWA,1,WA1,WA2,WA3)
C
C ON THE FIRST ITERATION AND IF MODE IS 1, SCALE ACCORDING
C TO THE NORMS OF THE COLUMNS OF THE INITIAL JACOBIAN.
C
IF (ITER .NE. 1) GO TO 70
IF (MODE .EQ. 2) GO TO 50
DO 40 J = 1, N
DIAG(J) = WA2(J)
IF (WA2(J) .EQ. ZERO) DIAG(J) = ONE
40 CONTINUE
50 CONTINUE
C
C ON THE FIRST ITERATION, CALCULATE THE NORM OF THE SCALED X
C AND INITIALIZE THE STEP BOUND DELTA.
C
DO 60 J = 1, N
WA3(J) = DIAG(J)*X(J)
60 CONTINUE
XNORM = ENORM(N,WA3)
DELTA = FACTOR*XNORM
IF (DELTA .EQ. ZERO) DELTA = FACTOR
70 CONTINUE
C
C FORM (Q TRANSPOSE)*FVEC AND STORE IN QTF.
C
DO 80 I = 1, N
QTF(I) = FVEC(I)
80 CONTINUE
DO 120 J = 1, N
IF (FJAC(J,J) .EQ. ZERO) GO TO 110
SUM = ZERO
DO 90 I = J, N
SUM = SUM + FJAC(I,J)*QTF(I)
90 CONTINUE
TEMP = -SUM/FJAC(J,J)
DO 100 I = J, N
QTF(I) = QTF(I) + FJAC(I,J)*TEMP
100 CONTINUE
110 CONTINUE
120 CONTINUE
C
C COPY THE TRIANGULAR FACTOR OF THE QR FACTORIZATION INTO R.
C
SING = .FALSE.
DO 150 J = 1, N
L = J
JM1 = J - 1
IF (JM1 .LT. 1) GO TO 140
DO 130 I = 1, JM1
R(L) = FJAC(I,J)
L = L + N - I
130 CONTINUE
140 CONTINUE
R(L) = WA1(J)
IF (WA1(J) .EQ. ZERO) SING = .TRUE.
150 CONTINUE
C
C ACCUMULATE THE ORTHOGONAL FACTOR IN FJAC.
C
CALL QFORM(N,N,FJAC,LDFJAC,WA1)
C
C RESCALE IF NECESSARY.
C
IF (MODE .EQ. 2) GO TO 170
DO 160 J = 1, N
DIAG(J) = DMAX1(DIAG(J),WA2(J))
160 CONTINUE
170 CONTINUE
C
C BEGINNING OF THE INNER LOOP.
C
180 CONTINUE
C
C IF REQUESTED, CALL FCN TO ENABLE PRINTING OF ITERATES.
C
IF (NPRINT .LE. 0) GO TO 190
IFLAG = 0
IF (MOD(ITER-1,NPRINT) .EQ. 0)
* CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
IF (IFLAG .LT. 0) GO TO 300
190 CONTINUE
C
C DETERMINE THE DIRECTION P.
C
CALL DOGLEG(N,R,LR,DIAG,QTF,DELTA,WA1,WA2,WA3)
C
C STORE THE DIRECTION P AND X + P. CALCULATE THE NORM OF P.
C
DO 200 J = 1, N
WA1(J) = -WA1(J)
WA2(J) = X(J) + WA1(J)
WA3(J) = DIAG(J)*WA1(J)
200 CONTINUE
PNORM = ENORM(N,WA3)
C
C ON THE FIRST ITERATION, ADJUST THE INITIAL STEP BOUND.
C
IF (ITER .EQ. 1) DELTA = DMIN1(DELTA,PNORM)
C
C EVALUATE THE FUNCTION AT X + P AND CALCULATE ITS NORM.
C
IFLAG = 1
CALL FCN(N,WA2,WA4,FJAC,LDFJAC,IFLAG)
NFEV = NFEV + 1
IF (IFLAG .LT. 0) GO TO 300
FNORM1 = ENORM(N,WA4)
C
C COMPUTE THE SCALED ACTUAL REDUCTION.
C
ACTRED = -ONE
IF (FNORM1 .LT. FNORM) ACTRED = ONE - (FNORM1/FNORM)**2
C
C COMPUTE THE SCALED PREDICTED REDUCTION.
C
L = 1
DO 220 I = 1, N
SUM = ZERO
DO 210 J = I, N
SUM = SUM + R(L)*WA1(J)
L = L + 1
210 CONTINUE
WA3(I) = QTF(I) + SUM
220 CONTINUE
TEMP = ENORM(N,WA3)
PRERED = ZERO
IF (TEMP .LT. FNORM) PRERED = ONE - (TEMP/FNORM)**2
C
C COMPUTE THE RATIO OF THE ACTUAL TO THE PREDICTED
C REDUCTION.
C
RATIO = ZERO
IF (PRERED .GT. ZERO) RATIO = ACTRED/PRERED
C
C UPDATE THE STEP BOUND.
C
IF (RATIO .GE. P1) GO TO 230
NCSUC = 0
NCFAIL = NCFAIL + 1
DELTA = P5*DELTA
GO TO 240
230 CONTINUE
NCFAIL = 0
NCSUC = NCSUC + 1
IF (RATIO .GE. P5 .OR. NCSUC .GT. 1)
* DELTA = DMAX1(DELTA,PNORM/P5)
IF (DABS(RATIO-ONE) .LE. P1) DELTA = PNORM/P5
240 CONTINUE
C
C TEST FOR SUCCESSFUL ITERATION.
C
IF (RATIO .LT. P0001) GO TO 260
C
C SUCCESSFUL ITERATION. UPDATE X, FVEC, AND THEIR NORMS.
C
DO 250 J = 1, N
X(J) = WA2(J)
WA2(J) = DIAG(J)*X(J)
FVEC(J) = WA4(J)
250 CONTINUE
XNORM = ENORM(N,WA2)
FNORM = FNORM1
ITER = ITER + 1
260 CONTINUE
C
C DETERMINE THE PROGRESS OF THE ITERATION.
C
NSLOW1 = NSLOW1 + 1
IF (ACTRED .GE. P001) NSLOW1 = 0
IF (JEVAL) NSLOW2 = NSLOW2 + 1
IF (ACTRED .GE. P1) NSLOW2 = 0
C
C TEST FOR CONVERGENCE.
C
IF (DELTA .LE. XTOL*XNORM .OR. FNORM .EQ. ZERO) INFO = 1
IF (INFO .NE. 0) GO TO 300
C
C TESTS FOR TERMINATION AND STRINGENT TOLERANCES.
C
IF (NFEV .GE. MAXFEV) INFO = 2
IF (P1*DMAX1(P1*DELTA,PNORM) .LE. EPSMCH*XNORM) INFO = 3
IF (NSLOW2 .EQ. 5) INFO = 4
IF (NSLOW1 .EQ. 10) INFO = 5
IF (INFO .NE. 0) GO TO 300
C
C CRITERION FOR RECALCULATING JACOBIAN.
C
IF (NCFAIL .EQ. 2) GO TO 290
C
C CALCULATE THE RANK ONE MODIFICATION TO THE JACOBIAN
C AND UPDATE QTF IF NECESSARY.
C
DO 280 J = 1, N
SUM = ZERO
DO 270 I = 1, N
SUM = SUM + FJAC(I,J)*WA4(I)
270 CONTINUE
WA2(J) = (SUM - WA3(J))/PNORM
WA1(J) = DIAG(J)*((DIAG(J)*WA1(J))/PNORM)
IF (RATIO .GE. P0001) QTF(J) = SUM
280 CONTINUE
C
C COMPUTE THE QR FACTORIZATION OF THE UPDATED JACOBIAN.
C
CALL R1UPDT(N,N,R,LR,WA1,WA2,WA3,SING)
CALL R1MPYQ(N,N,FJAC,LDFJAC,WA2,WA3)
CALL R1MPYQ(1,N,QTF,1,WA2,WA3)
C
C END OF THE INNER LOOP.
C
JEVAL = .FALSE.
GO TO 180
290 CONTINUE
C
C END OF THE OUTER LOOP.
C
GO TO 30
300 CONTINUE
C
C TERMINATION, EITHER NORMAL OR USER IMPOSED.
C
IF (IFLAG .LT. 0) INFO = IFLAG
IFLAG = 0
IF (NPRINT .GT. 0) CALL FCN(N,X,FVEC,FJAC,LDFJAC,IFLAG)
RETURN
C
C LAST CARD OF SUBROUTINE HYBRJ.
C
END